#packages: 
library(dplyr)
library(ggplot2)
library(DeclareDesign)
library(texreg)
library(cobalt)
library(WeightIt)
library(texreg)

#reading in data and initial cleaning ---------------
data2 <- read.csv('/Users/amygrauley/Desktop/Notre Dame/CN Elite Cues Experiment/EXP 2 DATA FINAL.csv')

#how many people did not pass the attention check
count_somewhat_disagree_2 <- table(data2$attention == "Somewhat disagree")[TRUE]
count_somewhat_disagree_2

#recoding to a binary variable where "1" is passing the attention check and "0" is not 
data2 <- data2 %>%
  mutate(attention = case_when(
    attention == "Strongly disagree" ~ 0, 
    attention == "Somewhat disagree" ~ 1, 
    attention == "Neither agree nor disagree" ~ 0, 
    attention == "Somewhat agree" ~ 0, 
    attention == "Strongly agree" ~ 0
  ))

#remove all who did not pass attention check from the data 
data2 <- data2[data2$attention != "0", ]

count_nopass_2 <- table(data2$attention == "1")[TRUE]
count_nopass_2

#if needed
#remove first three empty rows from data2
data2 <- data2[!rownames(data2) %in% c("NA", "NA.1"), , drop = FALSE]


#creating bins for CN score groups 
data2 <- data2 %>% 
  mutate(CN.Score2 = case_when(CN.Score == 0 ~ 1, #rejecters
                               CN.Score > 0 & CN.Score < 0.5 ~ 2, #skeptics
                               CN.Score > 0.49 & CN.Score < 0.75 ~ 3, #sympathizers
                               CN.Score > 0.74 ~ 4)) #adherents
#does the same thing as the block ID category 



#Creating a 7-point PID scale 
data2 <- data2 %>%
  mutate(pid = case_when(
    party_id == "Republican" & r_strength == "A strong Republican" ~ 7,
    party_id == "Republican" & r_strength == "Not a very strong Republican" ~ 6,
    party_id == "Independent" & i_strength == "Republican" ~ 5,
    party_id == "Independent" & i_strength == "Neither" ~ 4,
    party_id == "Independent" & i_strength == "Democratic" ~ 3,
    party_id == "Democrat" & d_strength == "Not a very strong Democrat" ~ 2,
    party_id == "Democrat" & d_strength == "A strong Democrat" ~ 1,
    TRUE ~ NA_integer_
  ))

#Coding ideology as numeric
data2 <- data2 %>% 
  mutate(ideology = case_when(ideology == "Very conservative" ~ 7, 
                               ideology == "Conservative" ~ 6, 
                               ideology == "Slightly conservative " ~ 5, 
                               ideology == "Moderate" ~ 4, 
                               ideology == "Slightly liberal" ~ 3,
                               ideology == "Liberal" ~ 2, 
                               ideology == "Very liberal" ~ 1, )) 

#coding religious importance as numeric 
data2 <- data2 %>% 
  mutate(rel_importance = case_when(
    rel_importance == "Religion is the most important thing in my life " ~ 4, 
    rel_importance == "Religion is one among many important things in my life " ~ 3, 
    rel_importance == "Religion is not as important as other things in my life " ~ 2, 
    rel_importance == "Religion is not important in my life" ~ 1
  ))



#making CN_MEASUREMENTs into numeric data2
data2 <- data2 %>%
  mutate(cn_measurement_1 = case_when(
    cn_measurement_1 == "Strongly disagree" ~ 1, 
    cn_measurement_1 == "Somewhat disagree" ~ 2, 
    cn_measurement_1 == "Neither agree nor disagree" ~ 3, 
    cn_measurement_1 == "Somewhat agree" ~ 4, 
    cn_measurement_1 == "Strongly agree" ~ 5
  ),
  cn_measurement_2 = case_when(
    cn_measurement_2 == "Strongly disagree" ~ 1, 
    cn_measurement_2 == "Somewhat disagree" ~ 2, 
    cn_measurement_2 == "Neither agree nor disagree" ~ 3, 
    cn_measurement_2 == "Somewhat agree" ~ 4, 
    cn_measurement_2 == "Strongly agree" ~ 5
  ),
  cn_measurement_3 = case_when(
    cn_measurement_3 == "Strongly disagree" ~ 1, 
    cn_measurement_3 == "Somewhat disagree" ~ 2, 
    cn_measurement_3 == "Neither agree nor disagree" ~ 3, 
    cn_measurement_3 == "Somewhat agree" ~ 4, 
    cn_measurement_3 == "Strongly agree" ~ 5
  ),
  cn_measurement_4 = case_when(
    cn_measurement_4 == "Strongly disagree" ~ 1, 
    cn_measurement_4 == "Somewhat disagree" ~ 2, 
    cn_measurement_4 == "Neither agree nor disagree" ~ 3, 
    cn_measurement_4 == "Somewhat agree" ~ 4, 
    cn_measurement_4 == "Strongly agree" ~ 5
  ),
  cn_measurement_5 = case_when(
    cn_measurement_5 == "Strongly disagree" ~ 1, 
    cn_measurement_5 == "Somewhat disagree" ~ 2, 
    cn_measurement_5 == "Neither agree nor disagree" ~ 3, 
    cn_measurement_5 == "Somewhat agree" ~ 4, 
    cn_measurement_5 == "Strongly agree" ~ 5
  ))


#outcome variables -------

#coding vote choice as numeric 
data2 <- data2 %>%
  mutate(vote_choice = case_when(
    vote_choice == "Extremely unlikely" ~ 5, 
    vote_choice == "Somewhat unlikely" ~ 4, 
    vote_choice == "Neither likely nor unlikely" ~ 3, 
    vote_choice == "Somewhat likely" ~ 2, 
    vote_choice == "Extremely likely" ~ 1))

#coding William's party as numeric: 
data2 <- data2 %>%
  mutate(williams_party = case_when(
    williams_party == "Republican" ~ 1, 
    williams_party == "Independent" ~ 4, 
    williams_party == "Democrat" ~ 7))

#coding social ties measurement of CN as numeric , 5= most CN
data2 <- data2 %>%
  mutate(friend_event_cn. = case_when(
    friend_event_cn. == "Yes" ~ 5, 
    friend_event_cn. == "Maybe yes" ~ 4, 
    friend_event_cn. == "Unsure/Don't Know" ~ 3, 
    friend_event_cn. == "Maybe no" ~ 2, 
    friend_event_cn. == "No" ~ 1))

#coding conversational identification measurement of CN as numeric , 5= most CN
data2 <- data2 %>%
  mutate(cn_we = case_when(
    cn_we == "Always" ~ 5, 
    cn_we == "Often" ~ 4, 
    cn_we == "Seldom" ~ 3, 
    cn_we == "Never" ~ 2, 
    cn_we == "Other:" ~ 1))

#coding CN measurement as numeric, 5= most CN
data2 <- data2 %>%
  mutate(cn_identification = case_when( 
    cn_identification == "Strongly disagree" ~ 1, 
    cn_identification == "Somewhat disagree" ~ 2, 
    cn_identification == "Neither agree nor disagree" ~ 3, 
    cn_identification == "Somewhat agree" ~ 4, 
    cn_identification == "Strongly agree" ~ 5))

#coding public identification with researcher contact as numeric, 5= most CN
data2 <- data2 %>%
  mutate(further_conversation = case_when( 
    further_conversation == "Yes" ~ 5, 
    further_conversation == "No, but I do identify with Christian Nationalism" ~ 3, 
    further_conversation == "No, because I do not identify with Christian Nationalism" ~ 1))


#coding for balance tests

#coding religion as protestant/evangelical vs not, binary 1-0
data2 <- data2 %>%
  mutate(religion_binary = case_when( 
    religion == "Catholic" ~ 0, 
    religion == "Mainline Protestant" ~ 1, 
    religion == "Other Protestant" ~ 1, 
    religion == "Evangelical Protestant" ~ 1, 
    religion == "None of the above" ~ 0, 
    religion == "Other:" ~ 0, 
    religion == "Jewish" ~ 0, 
    religion == "Other Christian" ~ 0, 
    religion == "Latter-Day Saint" ~ 0, 
    religion == "Other Christian non-religious" ~ 0))

#coding gender as male/nonmale, binary 1-0
data2 <- data2 %>%
  mutate(gender_binary = case_when( 
    gender == "Man" ~ 1, 
    gender == "Non-binary / third gender-" ~ 0, 
    gender == "Prefer not to say" ~ 0, 
    gender == "Woman" ~ 0))

data2 <- data2 %>%
  mutate(income_numeric = case_when( 
    income == "Less than $10,000" ~ 1, 
    income == "$10,000 - $19,999" ~ 2, 
    income == "$20,000 - $29,999" ~ 3, 
    income == "$30,000 - $39,999" ~ 4, 
    income ==  "$40,000 - $49,999" ~ 5, 
    income == "$50,000 - $59,999" ~ 6, 
    income == "$60,000 - $69,999" ~ 7, 
    income ==  "$70,000 - $79,999" ~ 8, 
    income ==  "$80,000 - $89,999" ~ 9, 
    income ==  "$90,000 - $99,999" ~ 10,
    income == "$100,000 - $149,999" ~ 11, 
    income == "More than $150,000" ~12))

data2 <- data2 %>%
  mutate(age_numeric = case_when( 
    age == "18 - 24" ~ 1, 
    age == "25 - 34" ~ 2, 
    age == "35 - 44" ~ 3, 
    age == "45 - 54" ~ 4, 
    age ==  "55 - 64" ~ 5, 
    age == "65 - 74" ~ 6, 
    age == "75 - 84" ~ 7))

#creating a binary variable for race
data2 <- data2 %>% 
  mutate(race_binary = case_when(
    race == "White, non- Hispanic"  ~ 1, 
    race == "Hispanic " ~ 0, 
    race == "Asian" ~ 0, 
    race == "Black or African American, non- Hispanic" ~ 0, 
    race == "American Indian or Alaskan Native" ~ 0, 
    race ==  "Multiple races, non-Hispanic" ~ 0, 
    race == "Other:" ~ 0, 
    race == "Native Hawaiian or other Pacific Islander" ~ 0
  ))


#turnCN.Score to numeric
data2$CN.Score <- as.numeric(data2$CN.Score)

#results for the main paper: no controls, mediated by CN.Score----------------

# Fit the model
model_nc_2 <- estimatr::lm_robust(cn_identification ~ Treatment + as.numeric(CN.Score),
                                data = data2)
summary(model_nc_2)

#LHP for two sets of treatments
linhyp_model_nc_LC <- car::linearHypothesis(model_nc_2, "TreatmentLCP = TreatmentLCN")
linhyp_model_nc_LC

linhyp_model_nc_RL <- car::linearHypothesis(model_nc_2, "TreatmentRLP = TreatmentRLN")
linhyp_model_nc_RL


#output tables
# Positive vs negative - Ha1 and Ha2
main_list_nc_2 <- list(model_nc_2)
texreg(main_list_nc_2, include.ci = FALSE)





#public ID-----------------------

#Creating scale for public identification: using questions FRIEND_EVENT_CN, CN_WE, FURTHER_CONVERSATION
data2$public_id_scale2 <-  rowMeans(data2[, c("friend_event_cn.", "cn_we", "further_conversation")], na.rm = TRUE)

#without friend_event_cn
data2$public_id_scale_2 <-  rowMeans(data2[, c("cn_we", "further_conversation")], na.rm = TRUE)

#public results without controls 

# Run the linear regression model
model_public_nc_2 <- estimatr::lm_robust(public_id_scale2 ~ Treatment+ as.numeric(CN.Score), 
                                       data = data2)
summary(model_public_nc_2)

linhyp_public_nc_RL <- car::linearHypothesis(model_public_nc_2, "TreatmentRLP = TreatmentRLN")
linhyp_public_nc_RL

linhyp_public_nc_LC <- car::linearHypothesis(model_public_nc_2, "TreatmentLCP = TreatmentLCN")
linhyp_public_nc_LC

#without friend metric in scale 
model_public_nc_2.1 <- estimatr::lm_robust(public_id_scale_2 ~ Treatment+ as.numeric(CN.Score), 
                                        data = data2)
summary(model_public_nc_2.1)


#output tables
# Positive vs negative - Ha1 and Ha2
public_list_nc_2 <- list(model_public_nc_2)
texreg(public_list_nc_2, include.ci = FALSE)


#all models in one list -----------
alllist <- list(model_nc, model_public_nc, model_nc_2, model_public_nc_2)
texreg(alllist, include.ci=FALSE)



#loading in data from durability check ---------------

data2_dur <- read.csv('/Users/amygrauley/Desktop/Notre Dame/CN Elite Cues Experiment/EXP2 DUR DATA FINAL.csv')

# Rename the column to match for merging 
colnames(data2_dur)[colnames(data2_dur) == "ID_DUR"] <- "id"

#coding social ties measurement of CN as numeric , 5= most CN
data2_dur <- data2_dur %>%
  mutate(DURABILITY_1. = case_when(
    DURABILITY_1. == "Yes" ~ 5, 
    DURABILITY_1. == "Maybe yes" ~ 4, 
    DURABILITY_1. == "Unsure/Don't Know" ~ 3, 
    DURABILITY_1. == "Maybe no" ~ 2, 
    DURABILITY_1. == "No" ~ 1))

#coding conversational identification measurement of CN as numeric , 5= most CN
data2_dur <- data2_dur %>%
  mutate(DURABILITY_2 = case_when(
    DURABILITY_2 == "Always" ~ 5, 
    DURABILITY_2 == "Often" ~ 4, 
    DURABILITY_2 == "Seldom" ~ 3, 
    DURABILITY_2 == "Never" ~ 2, 
    DURABILITY_2 == "Other:" ~ 1))

#coding CN measurement as numeric, 5= most CN
data2_dur <- data2_dur %>%
  mutate(DURABILITY_3 = case_when( 
    DURABILITY_3 == "Strongly disagree" ~ 1, 
    DURABILITY_3 == "Somewhat disagree" ~ 2, 
    DURABILITY_3 == "Neither agree nor disagree" ~ 3, 
    DURABILITY_3 == "Somewhat agree" ~ 4, 
    DURABILITY_3 == "Strongly agree" ~ 5))

#coding public identification with researcher contact as numeric, 5= most CN
data2_dur <- data2_dur %>%
  mutate(DURABILITY_4 = case_when( 
    DURABILITY_4 == "Yes" ~ 5, 
    DURABILITY_4 == "No, but I do identify with Christian Nationalism" ~ 3, 
    DURABILITY_4 == "No, because I do not identify with Christian Nationalism" ~ 1))

#adding the treatment into the data frame
#creating a data frame that pulls ID and treatment assignment from data_correct 
# Subset the data_correct data frame to select only the required columns
new_columns <- data2[, c("id", "Treatment", "CN.Score", "public_id_scale_2", "race", "pid",
                         "rel_importance", "gender", "ideology", "income", "age", "CN.Score", "religion", 
                         "CN.Score2")]

# Merge  with data2_dur based on participantId
data2_dur <- merge(data2_dur, new_columns, by = "id")

data2_dur <- na.omit(data2_dur)

#models for durability check --------------
# Fit the model

dur_check2 <- estimatr::lm_robust(DURABILITY_3 ~ Treatment + as.numeric(CN.Score),
                                 data = data2_dur)
summary(dur_check2)

linhyp_nc_RL_dur <- car::linearHypothesis(dur_check2, "TreatmentRLP = TreatmentRLN")
linhyp_nc_RL_dur

linhyp_nc_LC_dur <- car::linearHypothesis(dur_check2, "TreatmentLCP = TreatmentLCN")
linhyp_nc_LC_dur



#public ID durability check

dur_check2_pub <- estimatr::lm_robust(public_id_scale_2 ~ Treatment + as.numeric(CN.Score),
                                  data = data2_dur)
summary(dur_check2_pub)

linhyp_nc_RL_dur_pub <- car::linearHypothesis(dur_check2_pub, "TreatmentRLP = TreatmentRLN")
linhyp_nc_RL_dur_pub

linhyp_public_nc_LC_dur_pub <- car::linearHypothesis(dur_check2_pub, "TreatmentLCP = TreatmentLCN")
linhyp_public_nc_LC_dur_pub

#all durability models in one list 
alllist_dur <- list(dur_check, dur_check_pub, dur_check2,  dur_check2_pub)
texreg(alllist_dur, include.ci=FALSE)



#coefficient plots for both experiments --------------------
#active id
# Impute missing values with the mean for Experiment 1
mean_CN_IDENTIFICATION <- mean(data_correct$CN_IDENTIFICATION, na.rm = TRUE)
data_correct$CN_IDENTIFICATION[is.na(data_correct$CN_IDENTIFICATION)] <- mean_CN_IDENTIFICATION

# Impute missing values with the mean for Experiment 2
mean_cn_identification_exp2 <- mean(data2$cn_identification, na.rm = TRUE)
data2$cn_identification[is.na(data2$cn_identification)] <- mean_cn_identification_exp2

# Remove rows with NA values and ensure valid Treatments for plotting
plot_data_exp1 <- data_correct %>% filter(!is.na(CN_IDENTIFICATION) & !is.na(Treatment))
plot_data_exp2 <- data2 %>% filter(!is.na(cn_identification) & !is.na(Treatment) & Treatment != "")

# Normalize CN_IDENTIFICATION to a 0-1 scale for Experiment 1
plot_data_exp1 <- plot_data_exp1 %>%
  mutate(normalized_CN_IDENTIFICATION = (CN_IDENTIFICATION - min(CN_IDENTIFICATION)) / (max(CN_IDENTIFICATION) - min(CN_IDENTIFICATION)))

# Normalize cn_identification to a 0-1 scale for Experiment 2
plot_data_exp2 <- plot_data_exp2 %>%
  mutate(normalized_cn_identification = (cn_identification - min(cn_identification)) / (max(cn_identification) - min(cn_identification)))

# Combine data for facet plotting
summary_data_exp1 <- plot_data_exp1 %>%
  group_by(Treatment) %>%
  summarise(mean_value = mean(normalized_CN_IDENTIFICATION),
            lower_ci = mean_value - qt(0.975, length(normalized_CN_IDENTIFICATION)) * sd(normalized_CN_IDENTIFICATION) / sqrt(length(normalized_CN_IDENTIFICATION)),
            upper_ci = mean_value + qt(0.975, length(normalized_CN_IDENTIFICATION)) * sd(normalized_CN_IDENTIFICATION) / sqrt(length(normalized_CN_IDENTIFICATION)),
            Experiment = "Experiment 1")

summary_data_exp2 <- plot_data_exp2 %>%
  group_by(Treatment) %>%
  summarise(mean_value = mean(normalized_cn_identification),
            lower_ci = mean_value - qt(0.975, length(normalized_cn_identification)) * sd(normalized_cn_identification) / sqrt(length(normalized_cn_identification)),
            upper_ci = mean_value + qt(0.975, length(normalized_cn_identification)) * sd(normalized_cn_identification) / sqrt(length(normalized_cn_identification)),
            Experiment = "Experiment 2")

combined_summary_data <- bind_rows(summary_data_exp1, summary_data_exp2)

# Ensure the correct order of Treatment levels
combined_summary_data$Treatment <- factor(combined_summary_data$Treatment, 
                                          levels = c("P", "N", "C", "LCN", "LCP", "RLN", "RLP"))

# Create the plot
ggplot(combined_summary_data, aes(x = mean_value, y = Treatment)) +
  geom_point(position = position_dodge(width = 0.2), color = "black", size = 4, stroke = 0) +
  geom_errorbarh(aes(xmin = lower_ci, xmax = upper_ci), height = 0.3, color = "black") +
  labs(x = "Normalized Mean of CN Identification", y = "Treatment") +
  theme_minimal() +
  scale_y_discrete(labels = c("P" = "Positive Cue", "N" = "Negative Cue", 
                              "C" = "Control", "LCN" = "Love for Country- Negative", 
                              "LCP" = "Love for Country- Positive", 
                              "RLN" = "Religious Liberty- Negative", 
                              "RLP" = "Religious Liberty- Positive"), 
                   expand = c(0.2, 0)) +
  scale_x_continuous(limits = c(0, 0.5)) +
  geom_text(aes(label = round(mean_value, 3)),
            position = position_dodge(width = 0.2),
            vjust = 2,
            size = 3,
            color = "black") +
  facet_wrap(~Experiment, ncol = 2) +
  theme(panel.spacing = unit(2, "lines"))


# Create faceted dot and whisker plot with dividing line
ggplot(combined_summary_data, aes(x = mean_value, y = Treatment)) +
  geom_point(position = position_dodge(width = 0.2), color = "black", size = 4, stroke = 0) +
  geom_errorbarh(aes(xmin = lower_ci, xmax = upper_ci), height = 0.3, color = "black") +
  labs(x = "Normalized Mean of CN Identification", y = "Treatment") +
  theme_minimal() +
  scale_y_discrete(labels = c("P" = "Positive Cue (1)", "N" = "Negative Cue (1)", 
                              "C" = "Control", "LCN" = "Love for Country- Negative (2)", 
                              "LCP" = "Love for Country- Positive (2)", 
                              "RLN" = "Religious Liberty- Negative (2) ", 
                              "RLP" = "Religious Liberty- Positive (2) "), 
                   expand = c(0.2, 0)) +
  scale_x_continuous(limits = c(0, 0.5)) +
  geom_text(aes(label = round(mean_value, 3)),
            position = position_dodge(width = 0.2),
            vjust = 2,
            size = 3,
            color = "black") +
  facet_wrap(~Experiment, ncol = 2) +
  theme(panel.spacing = unit(2, "lines"))




#public id 
# Impute missing values with the mean for Experiment 1
mean_public_id_scale <- mean(data_correct$public_id_scale, na.rm = TRUE)
data_correct$public_id_scale[is.na(data_correct$public_id_scale)] <- mean_public_id_scale

# Impute missing values with the mean for Experiment 2
mean_public_id_exp2 <- mean(data2$public_id_scale2, na.rm = TRUE)
data2$public_id_scale2[is.na(data2$public_id_scale2)] <- mean_public_id_exp2

# Remove rows with NA values and ensure valid Treatments for plotting
plot_data_exp1_pub <- data_correct %>% filter(!is.na(public_id_scale) & !is.na(Treatment))
plot_data_exp2_pub <- data2 %>% filter(!is.na(public_id_scale2) & !is.na(Treatment) & Treatment != "")

# Normalize public_id_scale to a 0-1 scale for Experiment 1
plot_data_exp1_pub <- plot_data_exp1_pub %>%
  mutate(normalized_pub_id = (public_id_scale - min(public_id_scale)) / 
           (max(public_id_scale) - min(public_id_scale)))

# Normalize public_id_scale2 to a 0-1 scale for Experiment 2
plot_data_exp2_pub <- plot_data_exp2_pub %>%
  mutate(normalized_pub_id = (public_id_scale2 - min(public_id_scale2)) / 
           (max(public_id_scale2) - min(public_id_scale2)))

# Combine data for facet plotting
summary_data_exp1_pub <- plot_data_exp1_pub %>%
  group_by(Treatment) %>%
  summarise(
    mean_value = mean(normalized_pub_id),
    sd_value = sd(normalized_pub_id),
    n = n(),
    lower_ci = ifelse(n > 1,
                      mean_value - qt(0.975, df = n - 1) * sd_value / sqrt(n),
                      NA),
    upper_ci = ifelse(n > 1,
                      mean_value + qt(0.975, df = n - 1) * sd_value / sqrt(n),
                      NA),
    Experiment = "Experiment 1"
  )

summary_data_exp2_pub <- plot_data_exp2_pub %>%
  group_by(Treatment) %>%
  summarise(
    mean_value = mean(normalized_pub_id),
    sd_value = sd(normalized_pub_id),
    n = n(),
    lower_ci = ifelse(n > 1,
                      mean_value - qt(0.975, df = n - 1) * sd_value / sqrt(n),
                      NA),
    upper_ci = ifelse(n > 1,
                      mean_value + qt(0.975, df = n - 1) * sd_value / sqrt(n),
                      NA),
    Experiment = "Experiment 2"
  )



combined_summary_data_pub <- bind_rows(summary_data_exp1_pub, summary_data_exp2_pub)

# Ensure the correct order of Treatment levels
combined_summary_data_pub$Treatment <- factor(combined_summary_data_pub$Treatment, 
                                          levels = c("P", "N", "C", "LCN", "LCP", "RLN", "RLP"))

# Create the plot
ggplot(combined_summary_data_pub, aes(x = mean_value, y = Treatment)) +
  geom_point(position = position_dodge(width = 0.2), color = "black", size = 4, stroke = 0) +
  geom_errorbarh(aes(xmin = lower_ci, xmax = upper_ci), height = 0.3, color = "black") +
  labs(x = "Normalized Mean of Public CN Identification", y = "Treatment") +
  theme_minimal() +
  scale_y_discrete(labels = c("P" = "Positive Cue", "N" = "Negative Cue", 
                              "C" = "Control", "LCN" = "Love for Country- Negative", 
                              "LCP" = "Love for Country- Positive", 
                              "RLN" = "Religious Liberty- Negative", 
                              "RLP" = "Religious Liberty- Positive"), 
                   expand = c(0.2, 0)) +
  scale_x_continuous(limits = c(0, 0.5)) +
  geom_text(aes(label = round(mean_value, 3)),
            position = position_dodge(width = 0.2),
            vjust = 2,
            size = 3,
            color = "black") +
  facet_wrap(~Experiment, ncol = 2) +
  theme(panel.spacing = unit(2, "lines"))


# Create faceted dot and whisker plot with dividing line
ggplot(combined_summary_data_pub, aes(x = mean_value, y = Treatment)) +
  geom_point(position = position_dodge(width = 0.2), color = "black", size = 4, stroke = 0) +
  geom_errorbarh(aes(xmin = lower_ci, xmax = upper_ci), height = 0.3, color = "black") +
  labs(x = "Normalized Mean of PublicCN Identification", y = "Treatment") +
  theme_minimal() +
  scale_y_discrete(labels = c("P" = "Positive Cue (1)", "N" = "Negative Cue (1)", 
                              "C" = "Control", "LCN" = "Love for Country- Negative (2)", 
                              "LCP" = "Love for Country- Positive (2)", 
                              "RLN" = "Religious Liberty- Negative (2) ", 
                              "RLP" = "Religious Liberty- Positive (2) "), 
                   expand = c(0.2, 0)) +
  scale_x_continuous(limits = c(0, 0.52))+
  geom_text(aes(label = round(mean_value, 3)),
            position = position_dodge(width = 0.2),
            vjust = 2,
            size = 3,
            color = "black") +
  facet_wrap(~Experiment, ncol = 2) +
  theme(panel.spacing = unit(2, "lines"))

